//
// Created by Administrator on 2021/11/23.
//

#include "RayCasting.h"
#include "Utils.h"

RayCasting::RayCasting(Volume *volume) {
    this->volume = volume;
    rgba = new float[400 * 500 * 4];
    rotation_matrix();
}

void RayCasting::rotation_matrix() {
    // 图像空间的基向量
    float XX[3], YY[3], ZZ[3];
    ZZ[0] = eye[0] - center[0];
    ZZ[1] = eye[1] - center[1];
    ZZ[2] = eye[2] - center[2];
    Utils::cross_prod(XX, up, ZZ);
    Utils::cross_prod(YY, ZZ, XX);
    Utils::normalize(XX, XX);
    Utils::normalize(YY, YY);
    Utils::normalize(ZZ, ZZ);
    // 由图像空间基向量构成旋转矩阵
    R[0] = XX[0];
    R[1] = YY[0];
    R[2] = ZZ[0];
    R[3] = XX[1];
    R[4] = YY[1];
    R[5] = ZZ[1];
    R[6] = XX[2];
    R[7] = YY[2];
    R[8] = ZZ[2];
}

void RayCasting::composite(float *temp, int x, int y) {
    // 采样步长
    int stepSize = 1;
    // 累计颜色值
    float cumColor[4]{0.0, 0.0, 0.0, 0.0};
    // 投射光线起点、方向
    float pos[3], dir[3];
    // 光线与包围盒近视点处的交点坐标
    float startPos[3];
    // 采样点坐标
    float samplePos[3];
    // 采样点颜色
    float sampleColor[4];
    // 采用平行投影，故在图像空间中投射光线的方向(0,0,-1),起点(x,y,0)
    pos[0] = x;
    pos[1] = y;
    pos[2] = 0;

    // 将光线描述转换到物体空间
    // 光线方向在物体空间的表达
    dir[0] = -R[2];
    dir[1] = -R[5];
    dir[2] = -R[8];
    // 旋转
    Utils::matrix_mul_vec(pos, R, pos);
    // 平移
    pos[0] += T[0];
    pos[1] += T[1];
    pos[2] += T[2];

    // 判断光线与包围盒是否相交
    if (intersection(startPos, pos, dir)) {
        samplePos[0] = startPos[0];
        samplePos[1] = startPos[1];
        samplePos[2] = startPos[2];
        // 当光线射出包围盒或累计不透明度超过 1 时中止合成
        while (check_in_box(samplePos) && cumColor[3] < 1) {
            // 三线性插值获得采样点处的颜色及不透明度
            tr_interpolation(sampleColor, samplePos);
            // 合成颜色及不透明度,采用的是从前到后的合成公式
            cumColor[0] += sampleColor[0] * sampleColor[3] * (1 - cumColor[3]);  // R
            cumColor[1] += sampleColor[1] * sampleColor[3] * (1 - cumColor[3]);  // G
            cumColor[2] += sampleColor[2] * sampleColor[3] * (1 - cumColor[3]);  // B
            cumColor[3] += sampleColor[3] * (1 - cumColor[3]);                   // A
            // 下一个采样点
            samplePos[0] += dir[0] * stepSize;
            samplePos[1] += dir[1] * stepSize;
            samplePos[2] += dir[2] * stepSize;
        }
        temp[0] = cumColor[0];
        temp[1] = cumColor[1];
        temp[2] = cumColor[2];
        temp[3] = cumColor[3];
        return;
    }
    temp[0] = temp[1] = temp[2] = temp[3] = 1.0;  //若光线与包围盒不相交，赋白色
}

bool RayCasting::check_in_box(const float *point) {
    auto dim = volume->get_dim();

    if (point[0] < 0 || point[0] >= dim[0] || point[1] < 0 ||
        point[1] >= dim[1] || point[2] < 0 || point[2] >= dim[2])
        return false;
    else
        return true;
}

bool RayCasting::intersection(float *startPos, const float *pos, const float *dir) {
    auto dim = volume->get_dim();
    float nearScale = -1000000;
    float scale1, scale2;
    // 光线与包围盒平行于YOZ的2个平面交点
    if ((dir[0] <= -EPSILON) || (dir[0] >= EPSILON)) {
        scale1 = (0 - pos[0]) / dir[0];
        scale2 = (dim[0] - 1 - pos[0]) / dir[0];
        //选出靠近视点的交点，并与当前候选点比较，保留较远者
        if (scale1 < scale2) {
            if (scale1 > nearScale) nearScale = scale1;
        } else {
            if (scale2 > nearScale) nearScale = scale2;
        }
    }
    // 光线与包围盒平行于ZOX的2个平面交点
    if ((dir[1] <= -EPSILON) || (dir[1] >= EPSILON)) {
        scale1 = (0 - pos[1]) / dir[1];
        scale2 = (dim[1] - 1 - pos[1]) / dir[1];
        //选出靠近视点的交点，并与当前候选点比较，保留较远者
        if (scale1 < scale2) {
            if (scale1 > nearScale) nearScale = scale1;
        } else {
            if (scale2 > nearScale) nearScale = scale2;
        }
    }
    // 光线与包围盒平行于XOY的2个平面交点
    if ((dir[2] <= -EPSILON) || (dir[2] >= EPSILON)) {
        scale1 = (0 - pos[2]) / dir[2];
        scale2 = (dim[2] - 1 - pos[2]) / dir[2];
        //选出靠近视点的交点，并与当前候选点比较，保留较远者
        if (scale1 < scale2) {
            if (scale1 > nearScale) nearScale = scale1;
        } else {
            if (scale2 > nearScale) nearScale = scale2;
        }
    }

    startPos[0] = pos[0] + nearScale * dir[0];
    startPos[1] = pos[1] + nearScale * dir[1];
    startPos[2] = pos[2] + nearScale * dir[2];

    // 判断该点是否在包围盒内
    return check_in_box(startPos);
}

void RayCasting::tr_interpolation(float *sampleColor, const float *pos) {
    auto dim = volume->get_dim();
    auto cData = volume->get_c_data();

    int x0, y0, z0, x1, y1, z1;
    float fx, fy, fz;
    float v0, v1, v2, v3, v4, v5, v6;
    int Slicesize = dim[0] * dim[1] * 4;
    int Stepsize = dim[0] * 4;
    x0 = (int) pos[0];  //整数部分
    y0 = (int) pos[1];
    z0 = (int) pos[2];
    fx = pos[0] - x0;  //小数部分
    fy = pos[1] - y0;
    fz = pos[2] - z0;
    x1 = x0 + 1;
    y1 = y0 + 1;
    z1 = z0 + 1;
    // 防止越界
    if (x1 >= dim[0]) x1 = dim[0] - 1;
    if (y1 >= dim[1]) y1 = dim[1] - 1;
    if (z1 >= dim[2]) z1 = dim[2] - 1;

    for (int i = 0; i < 4; i++) {
        // 采样点处的值由邻近的8个点插值获得
        v0 = cData[z0 * Slicesize + y0 * Stepsize + 4 * x0 + i] * (1 - fx) +
             cData[z0 * Slicesize + y0 * Stepsize + 4 * x1 + i] * fx;
        v1 = cData[z0 * Slicesize + y1 * Stepsize + 4 * x0 + i] * (1 - fx) +
             cData[z0 * Slicesize + y1 * Stepsize + 4 * x1 + i] * fx;
        v2 = cData[z1 * Slicesize + y0 * Stepsize + 4 * x0 + i] * (1 - fx) +
             cData[z1 * Slicesize + y0 * Stepsize + 4 * x1 + i] * fx;
        v3 = cData[z1 * Slicesize + y1 * Stepsize + 4 * x0 + i] * (1 - fx) +
             cData[z1 * Slicesize + y1 * Stepsize + 4 * x1 + i] * fx;
        v4 = v0 * (1 - fy) + v1 * fy;
        v5 = v2 * (1 - fy) + v3 * fy;
        v6 = v4 * (1 - fz) + v5 * fz;
        // 防止越界
        if (v6 > 1) v6 = 1;
        sampleColor[i] = v6;
    }
}

RayCasting::~RayCasting() {
    delete[]rgba;
}
